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A dynamical theory of geophysical precipitation pattern formation is presented and applied to 
irreversible calcium carbonate (travertine) deposition. Specific systems studied here are the terraces 
and domes observed at geothermal hot springs, such as those at Yellowstone National Park, and 
speleothems, particularly stalactites and stalagmites. The theory couples the precipitation front 
dynamics with shallow water flow, including corrections for turbulent drag and curvature effects. 
In the absence of capillarity and with a laminar flow proflle, the theory predicts a one-parameter 
family of steady state solutions to the moving boundary problem describing the precipitation front. 
These shapes match well the measured shapes near the vent at the top of observed travertine domes. 
Closer to the base of the dome, the solutions deviate from observations, and circular symmetry is 
broken by a fluting pattern, which we show is associated with capillary forces causing thin fllm 
break-up. We relate our model to that recently proposed for stalactite growth, and calculate the 
linear stability spectrum of both travertine domes and stalactites. Lastly, we apply the theory to the 
problem of precipitation pattern formation arising from turbulent flow down an inclined plane, and 
identify a linear instability that underlies scale-invariant travertine terrace formation at geothermal 
hot springs. 

PACS numbers: 05.45.Ra, 87.23.n, 47.54.-r, 89.75.Kd, 47.20. Hw, 47.15.gm, 47.55.np 



I. INTRODUCTION 

Geophysical pattern formation concerns how geolog- 
ical patterns and landscapes are formed as a result of 
the underlying physical and chemical dynamics. The 
aim is to predict the static, dynamical and statisti- 
cal properties of the variety of geological structures 
formed. Recently studied examples include travertine 
motifs, namely damspj, domes [2] and terraces [21 El HI [5], 
stalactites [SI |7], as well as that of other patterns such as 
sand dunes [HI H], black smoker chimneys at hydrother- 
mal vents [To] . columnar joints [TT] and braided river 
networks [12]. 

This paper focuses on the formation of travertine struc- 
tures near geothermal hot springs. In such systems, hot 
spring water emerges from a vent, and deposits calcium 
carbonate as a mineral generally termed travertine as it 
degasses carbon dioxide [H [31 H]. The formation of 
stalactites in limestone caves, which are also caused by 
carbonate precipitation, will also be briefly discussed. 

The majority of the work done on the subject has fo- 
cused on the microscopic aspects of the problem, such 
as the role of biomineralization due to thermophilic 
microbes[Sl H], the CO2 degassing mechanisms [Ol [Hj, 
mineral compositions [151 US] and crystal structure [TT] 
[T5] . Here we are interested in the formation of macro- 
scopic structures and motifs, such as domes, stalactites, 
and terraces[21, which are universal, i.e., independent 
of microscopic details. In addition, we are interested 
in the resulting patterns and their correlations, rather 
than absolute rates of growth; accordingly, microscopic 
mechanisms that contribute to kinetics, including nucle- 
ation processes and potential biomineralization effects, 
are present in our work through the choice of time scale. 




FIG. 1: (Color online) Travertine formation at Angel Terrace, 
Mammoth Hot Springs, WY, showing a large pond, of order 
1 meter in diameter, and smaller features. 

There are no extra terms in the equations of motion 
whose presence can be attributed specially to any one 
of these microscopic processes. 

There are two principal mathematical difficulties en- 
countered in studying these macroscopic structures. 
First, the problem is highly nonlinear. As the carbon- 
ate is precipitated onto the surface, the surface evolves, 
which then changes the flow path of the fluid, thus affect- 
ing how precipitation takes place. This interplay between 
fluid flow and surface growth leads to a moving-boundary 
problem, which is mathematically difficult to solve. Sec- 
ond, the problem involves a variety of depositional pro- 
cesses, including solute advection, a complex sequence of 
chemical reactions, CO2 degassing, as well as mass trans- 
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fer between a solid and a liquid. Given that each of these 
processes is complicated and non-trivial to model on its 
own, a holistic approach capturing all of them would not 
be mathematically tractable. 

The purpose of this paper is to explore a simplified 
mathematical formulation of this problem that captures 
the essential large-scale dynamics. Because of the com- 
plexity of the problem, the resulting equations are very 
complicated, making it difficult, if not impossible, to un- 
derstand the whole flow system using this approach. It 
turns out, however, that the equations can be solved an- 
alytically under some simple situations, where symmetry 
can be exploited and simplifications can be made. The 
formations of domes [2^ and stalactites [BJ [7] are examples 
of such situations, as is the pioneering work of Wooding 
on travertine dams[T]. In these systems, there is a thin 
film of fiuid fiowing over the motif in a laminar fashion (in 
the case of domes and stalactites, for example). We will 
see that these simple motifs are straightforward to cal- 
culate in the case that capillary forces can be neglected. 
If the fluid fllm becomes too thin, due to its spreading 
over the surface, contact lines can be formed, resulting 
in rivulets and the breaking of pure rotational symme- 
try. In the case of domes, this is manifested in a fluting 
pattern near the base of the dome|lj. Such effects are 
difficult to include analytically, although we have previ- 
ously shown that they can be captured correctly using a 
cell dynamical system model [2], and this is discussed in 
more detail below. 

Although we cannot use this analytical theory to study 
the detailed shapes of the complex landscape of ponds 
and terraces, we are able to expose the dynamical linear 
instabilities, whose evolution into the nonlinear regime 
give rise to the landscape. We will see that the linear 
stability spectrum, in the absence of capillarity effects, 
always predicts a positive growth rate. The absence of a 
length scale arising in this calculation suggests that the 
actual landscapes might be scale invariant, a conclusion 
that is reinforced by our studies of the statistical proper- 
ties of these landscapes using our cell dynamical system 
model and photographic evidence p, [T^ 

The study reported here is a complement to our sim- 
ulation work [2] [ini [20] implemented as a cell dynamical 
system. This model has been shown to be capable of 
describing the actual dynamics |2], not only in the simple 
cases where the analytical approach is successful, but also 
in the fully nonlinear regime. For example, it has been 
shown that this cellular model generically gives rise to a 
complex, terraced landscape, which is similar to the one 
observed in the field. The cellular model also makes de- 
tailed predictions for the landscape statistics, including 
the pond area distribution and the distribution of pond 
anisotropy. In addition, the model successfully predicts 
that the main mode of pond or terrace growth is uphill 
pond inundation, a result confirmed by time-lapse pho- 
tographic studies. 

Although seemingly different, both the analytical ap- 
proach and the cell dynamical system approach incorpo- 



rate the same physics, and so should be expected to yield 
identical predictions. In [2] this was tested, by using the 
cellular model to solve the problem of dome formation. 
The analytical theory in the absence of surface tension 
cannot account for the fluting seen away from the vent of 
domes, because the fluting arises from contact line forma- 
tion. The analytical theory for domes, as we will discuss 
in detail below, contains one parameter that sets the scale 
for the patterns: this scale factor rg is a combination of 
the upward growth velocity, the mass transfer coefficient 
describing how material is incorporated into the growing 
substrate, the flux of water emerging from the vent, the 
gravitational acceleration and the fluid viscosity. When 
surface tension effects are included, the capillary length 
c?o must also be included. Thus, our theory is a two 
parameter theory for the entire range of travertine depo- 
sitional phenomena. The analytical theory can be used 
to predict the position on the dome at which capillary 
effects become important: this must occur at a location 
that is independent of the ratio ro/do, and hence this 
critical angle has a prescribed dependence on the un- 
derlying parameters which enter into the formula for tq. 
This prediction, arising from the analytical theory, was 
verifled to occur also in the computer simulations of the 
cellular model[2]. As a result, we conclude that the two 
formulations are indeed equivalent, and may be used in- 
terchangeably depending on which is more suited to the 
problem at hand. 

This paper is organized as follows. In Section [ll] 
we derive the equations governing the dynamics of fluid 
flow coupled to the moving boundary problem describing 
travertine precipitation. Section |III| describes the circu- 
larly symmetric solutions of these equations, and presents 
the linear stability analysis of the steady state uniformly 
translating solutions. We compare our analysis to a sim- 
ilar one[nj [7] that describes the shapes of stalactites in 
Section |IV| and compute the linear stability spectrum of 
these structures too. We turn in Section|Vl]to a study of 
turbulent flow down an inclined plane, and calculate the 
linear stability spectrum for the coupled flow and moving 
boundary problem, exposing the linear instability that is 
at the heart of the terraced landscape architecture. We 
conclude in Section fVIII 



II. MODEL FOR PRECIPITATION PATTERN 
FORMATION 

We consider a stream of water flowing over a terrain, 
from which calcium carbonate is then, due to geochemical 
processes to be discussed below, precipitated onto the 
landscape. The landscape is thus constantly changing 
in response to the fluid flow. This change of landscape, 
in turn, affects the flow path of the fluid, which than 
influences how subsequent precipitation takes place. We 
derive the governing equations describing both fluid flow 
and surface growth. We flrst focus on the surface growth, 
and related precipitation dynamics, and then move onto 
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the fluid flow. These two aspects will be combined to 
provide the complete description of the system. 



A. Surface Growth 

A surface can generally be characterized by the local 
curvature, k. In one dimension, or in cases where symme- 
try reduces the system to be effectively one dimensional, 
K is defined by 



K = 



de 



(1) 



where 9 is an angle between the local tangent of the curve 
and a fixed axis, and s is the arc length measured from 
some fixed point on the curve, as shown in Fig. ([2|. If 
the normal velocity u„ of the surface is prescribed ev- 
erywhere, then the evolution of the curvature follows the 
kinematic eauation |21l f22l [23] : 



Ok 
dt 
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(2) 



The time derivative in the equation is defined with re- 
spect to fixed 9. The first term in Eq. ^ describes the 
change in k due to the change in the overall scale of the 
object, whereas the second term describes the change in 
K at a point due to the difference in growth rates in the 
neighborhood of that point. 

Eq. ([2| is purely geometrical; for any given function 
Vn, the evolution of k is determined. So, physics enters 
in constructing a realistic and mathematically tractable 
model for ?;„, which, in the case considered here, de- 
pends on water chemistry, surface kinematics, chemical 
advection and fluid flow state. In carbonate systems, in 
additional to the CaCOa concentration, precipitation is 
mainly controlled by the CO2 concentration (partially 
reflected in the measured pH), which is also influenced 
by its temperature-dependent solubility in the fluid. As 
the pH increases or the temperature decreases, the solv- 
ability of CaCOs decreases and supersaturated CaCOa 
will be precipitated onto the surface. While the decrease 
in temperature is mainly due to heat loss to the envi- 
ronment, the increase in pH is due to the loss of CO2 
by a variety of outgassing mechanisms f 1 3| [T4]. Although 
the detailed water chemistry and depositional processes 
are quite complicated, for the purposes of the present 
work, it suffices to use a simplification of the governing 
chemical reactions: Gd?+ + 2HCO3 ^ CaC03(s) + H2O 
+ C02(g). In summary, the system tends to produce 
more CaCOs as CO2 concentration decreases through 
outgassing. 

Mass transfer between a fluid and a solid is a compli- 
cated problem,!, ,24, i25^; these nontrivial chemical pro- 
cesses only make it harder. A complete description of 
the precipitation dynamics, which will give us the nor- 
mal growth velocity Vn, involves writing down, in addi- 
tion to the fluid dynamics equations, advection-reaction- 
diffusion equations for each chemical and appropriate 



boundary conditions. Short et aZ. [SJIT] followed this ap- 
proach in the study of stalactite formation. What they 
found, after solving all these equations and taking limits 
appropriate for the timescales of interest to them, is that 
Vn is proportional to the local fluid thickness, h, with 
all the chemistry entering only into the proportionality 
constant. 

A simple interpretation of this result can be obtained 
by studying the scales of processes involved in stalac- 
tite formation, using parameter values from Ref. [7]. 
The fluid flow is a laminar flow, with Reynold's number 
about 0.01 — 1 . The thickness of the flow, h, is typically 
on the order of 10/im. The time scale for the concen- 
tration of CaCOa to equilibrate across the layer is thus 
h^/D ~ 0.1s, where D is the diffusion constant. Next, 
the traversal time, the time for a parcel of fluid to flow 
along the stalactites, is about 100s. Because only 1 per- 
cent of the total CaCOa mass is precipitated throughout 
the flow, we can assume that the CaCOa concentration, 
and thus the pH, are uniform both across the fluid layer 
and along the stalactite. The temperature can also be 
assumed to be constant since the fluid is so thin. The 
precipitation rate is then controlled only by the CaCOa 
available, which is proportional to the thickness of the 
fluid. 

In other carbonate systems, such as at travertine- 
forming hot springs, this relation between w„ and h does 
not hold simply due to the fact that the fluid thickness is 
larger, and the velocity is larger; as a result a turbulent 
boundary layer is formed near the precipitation front. 
What happens outside the boundary layer is too distant 
to affect precipitation near the boundary. In a turbu- 
lent flow, instead of depending on h, the precipitation 
front velocity depends on the fluid velocity [5S]. 
Wooding [1, in the study of steady-state dam formation, 
took this into consideration and arrived at the conclusion 
that Vn is directly proportional to the depth-averaged 
tangential fluid velocity, U, i.e. 



= GC/, 



(3) 



where G is a mass transfer coefficient depending on 
water chemistry and spectral features of the turbulent 
flow 1211 [13. For present purposes, the functional form 
of G is not of interest: we shall treat it as a phenomeno- 
logical parameter, and as we shall see, its role in the the- 
ory developed here is to contribute to the characteristic 
length scale tq of patterns. 

To summarize: all the details of water chemistry, 
including supersaturation, outgassing, solute diffusion, 
fluid turbulence, temperature and pH, which on their own 
are complicated processes and are nontrivial to model, 
enter into the picture only through a mass transfer co- 
efficient, G. In principle, G may exhibit spatial fluctua- 
tions; however, we shall assume that these are on a scale 
small compared to the features we are describing, and 
thus we will consider G to be a constant locally along 
the flow path. Over the entire geothermal spring system, 
it is possible that there will be a small spatial variation 
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FIG. 2: The coordinate system for the model of fluid flow 
coupled to precipitation moving boundary dynamics. 



where C/ is the Chezy coefhdent|26|, which empirically 
describes the energy lost due to turbulence, in a manner 
consistent with Kolmogorov's 1941 scaling theory of tur- 
bulence (K41)[251 [53], and s is the arc length measure 
from a reference point at the top, as shown in Fig. 

The de Saint- Venant equation only holds on flat sur- 
faces. When the surface grows, flow instabilities trigger 
various patterns to form; and the de Saint- Venant equa- 
tion is no longer vahd. For a general curved surface, the 
Dressier equation [30} [31] has to be used: 



1 duQ 
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where 



u(s, n, t) 
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in the mean value of G, but the weak dependence of G 
on governing parameters [T] |2H US] strongly suggests that 
this can reasonably be neglected. 



B. Fluid Dynamics 

A complete description of incompressible fluid dynam- 
ics is given by the Navier-Stokes equation 



E{s,t) C + hcose 



du 

m 
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u-Vu= — VP -I- lyV^u + g, 



(4) 



with V • M = for incompressibility, no-slip and stress- 
free boundary conditions at the solid-liquid and liquid- 
gas interfaces, respectively, where u, p, P, v and g arc 
the fluid velocity, density, pressure, viscosity and gravi- 
tational acceleration. We will use the Poiseuille solutions 
of the Navier-Stokes equations for domes, where the flow 
is laminar, but for turbulent flows, such as those which 
form the travertine terraces, we will employ a depth- 
averaging approximation, in conjunction with the Chezy 
approximation [26] for hydraulic friction. 

Since the spatial scale over which the landscape 
changes is usually much larger than the fluid thickness, 
i.e. hn <^ 1, we can make use of the shallow water ap- 
proximation and expand Eq. ^ in powers of /ik. If 
we take k to be zero, we arrive at the de Saint- Venant 
equation [57] 



d(Uh) dOJ'^h) Bh , . , 

with equation of continuity 

d{Uh) ^ 



dh 
'di 



ds 
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where C is the height of the underlying surface measured 
from a fixed horizontal axis, as shown in Fig. ([2|, ph 
is the pressure head at the fluid surface, p is the fluid 
density, E is the energy density and q is the local flux. 
When K is set to zero and is small, the Dressier equa- 
tions reduce to those of de Saint- Venant. 

As we have seen, the way fluid flows depends on the 
landscape it is flowing over, which itself is evolving over 
time. Now, Eq. (or Eq. Q) and Eq. ^ de- 

scribe these two dynamics, respectively. However, we do 
not have to consider both dynamics on the same footing 
because there is a separation of time scales; the rate of 
fluid flow is on the order of cm/sec, but the rate of precip- 
itation is on much slower geological scales. The latter is 
on the order of 1 mm/day and 1 cm/century in the cases 
of Yellowstone travertines [31 |3H [33] and stalactites [7], 
respectively. Accordingly the fluid flow responds quickly 
to the landscape, but the landscape responds extremely 
slowly to the fluid flow. We can then assume that the 
fluid flow is in its steady state when we discuss the land- 
scape evolution; i.e., we can drop all the time derivatives 
in the fluid flow equations. This quasi-stationary model 
will now be used to study the steady states of a variety 
of geological motifs and their stabihties. 



III. TRAVERTINE DOMES 
A. Steady state 

Our first example is the circularly symmetrical domes 
found in Yellowstone National Park, as shown in Fig. 



([3^). A number of approximations and simplifications 
can be made before we proceed. First, the growth rate 
of these domes is on the order of 1 — 5mm/day and the 
fluid flow rate is on the order of Imm/s, so we have a 
separation of time scales. Second, our field observations 
indicate that the thickness of the fluid film fiowing over 
the domes is very small compared to the curvature of 
the surface; thus, we make the approximation that the 
fluid is flowing down a (locally) constant slope. Third, 
as suggested by the field observations, the domes have a 
high degree of circular symmetry, so we can assume the 
solution to be circularly symmetrical and focus only on 
the radial part of the solution, which is effectively one 
dimensional. Fourth, the flow is apparently laminar, so 
we can use the Poiseuille-Hagen profile for the velocity 
in thin film: 



gh'^ sin 9 
2^ 



(12) 



where 9 is the slope of the surface and y is the transverse 
coordinate, as shown in Fig. (|2|. By assuming circular 
symmetry, h can be related to the axial distance from 
the vent, r, by the conservation of fluid volume: 



Q = l-nr j u{y)dy - 



sin 9 



(13) 



where Q is the total volumetric flux coming out of the 
vent. Eq. (12) and (13 1 can be combined to give 



U - 



u{y)dy 



a sm ( 



1/3 



(14) 



where a = gQ^ / \2'kv . We will see later that the assump- 
tion of laminar flow is self-consistently verified. Putting 
Eq. ([ll]) into Eq. ^ and using Eq. (^, gives 



dn 
dt 



G 
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1/3 



(15) 



This is the governing equation for the dome profile. Sug- 
gested by the shape of the dome, we seek a solution which 
steadily translates upwards without a change of shape, 
i.e., dtn\e = 0, with velocity vt. Eq. (15l gives 



G 



asin^ 



1/3 



Vt COS ( 



(16) 



Rearranging terms gives the shape of the dome as a one- 
parameter family of curves 



sin 9 



COS'' 



(17) 



where the scale factor ro = ^/CPaJvf. Eq. (17 1 is 
plotted in Fig. (|3b). Good agreement is obtained be- 



tween our theory and the observations below a critical 
angle 9c. From the fit, and the typical parameter values 
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FIG. 3: (Color online) Travertine dome at Mammoth Hot 
Springs, WY. (a) Dome whose central pond is 50.3cm in di- 
ameter, (b) Dome profile compared with theory and simula- 
tion of Ref. 2 . The black cmve is the analytical prediction 
from Eq. ( 17l, using the value ro — 43cm. The red filled cir- 



cles show the profile of a simulated dome, including the effects 
of surface tension. The blue dashed line is a consensus dome 
profile generated by averaging the dome shown with one other 
field observation. The blue filled squares show the profile of 
a simulated dome without surface tension|2j. 



G ^ 10~^, Vt ~ Imm/day and Q ~ Icm'^/sec, we obtain 
U ~ 25mm/sec and /i ~ 1 — 10mm, and a Reynolds num- 
ber. Re = Uh/v ^ 10 — 100. The assumption of laminar 
flow is self-consistently verified. 

The agreement between this analysis and observation 
shows that the growth of the dome is mainly determined 
by the geometry, because the only r dependence enters 
through the mass conservation, which is determined by 
geometry. To see this, suppose that the dome was a one 
dimensional object. Then, the mass conservation equa- 
tion, Eq. ( 13 ), would become Uh = qo, for some constant 
fiux qq, without any r dependence. Under the same ap- 



proximation of local fiatness, the final equation for U , Eq. 
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( 14 1, would thus be independent of r. We would then not 
be able to solve for r by substituting U into Eq. In 
this case, we would have to solve the equations without 
using the locally-flat approximation. In other words, the 
fact that we can ignore the details of the flow, by as- 
suming local flatness, to obtain the shape of the domes 
implies that geometry plays a more important role than 
fluid flow in the formation of domes. 

For angles 9 > 9c, the analytical profile deviates from 
our field photograph. The point of deviation is associ- 
ated with an apparent change in the dome morphology, 
with a fluting pattern superimposed on the dome pro- 
file. This is due to the effects of surface tension at the 
air- water-travertine interface [2] . Instead of covering the 
whole surface uniformly, the fiuid separates and covers 
only a fraction of the surface. Along the wetted surface, 
the regular growth law still applies and thus the surface 
grows, until a point at which the difference in heights 
between the wetted and dry surfaces is so large that the 
flow changes its path to flow along the dry surface. This 
process repeats itself and, on average, results in a slower 
growth when compared with a uniformly-wetted dome, 
so the theoretical prediction should be higher than the 
observation for 9 > 9c, as seen in Fig. ^p). The analyti- 
cal solution for the dome profile neglects surface tension, 
but leads to a prediction for the scaling dependence of 
the critical angle on the model parameters [2 . 

It is not trivial to include surface tension in our ana- 
lytical model, but its effect can be examined by using the 
cellular model, in which one can switch on and off surface 
tension. Fig. ^p), reproduced from Ref. |1] shows the 
prediction of dome shapes from the cellular model with 
and without surface tension. It is clear that by appropri- 
ate choice of do the simulation result coincides with the 
observation when surface tension is present, and agrees 
with the analytical prediction otherwise. This is direct 
evidence for the effect of surface tension near fluting. 

For completeness, we mention that this is not an ar- 
tifact of having "enough fitting parameters to fit an ele- 
phant" . In Ref. [2] was presented a scaling argument for 
the critical angle at which capillary effects become im- 
portant. The inclusion of surface tension introduces an 
additional length scale, namely, the capillary length, dc, 
into the problem. Now, the only other length scale in the 
problem is tq = \/ gG^Q"^ /i^vf . Since 9c is dimensionless, 
it can only depend on the ratio ro/dc and G. For a given 
chemical environment, G is fixed and we are left with the 
prediction, derived from our analytical solution, that 



= / 



(18) 



where f{x) is a scaling function. This data collapse, 
which predicts 9 depends not on the parameters sepa- 
rately, but only on the combination {gQ'^ / vvf) j dc, was 
verified using our discrete cellular model [2 , wherein the 
form of ]{x) was calculated. 



B. Linear stability analysis 

To complete the analysis, we study the stability of the 
solution, Eq. (17 1. By following the approach Liu and 



Goldenfeld used in studying the linear stability of den- 
dritic solidification [33], we consider a perturbed solution, 
riff) — f{9) + Sr{9)e^*, where r{9) is the solution in Eq. 
( 17 1 and Sr is a perturbation. Substituting this into the 



governing equation, Eq. (15), and expanding in Sr, we 
obtain 
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where the boundary conditions are 



5r{Q) = 0, 

for symmetric modes and 
d{5r{Q)) 



= 0, (19) 



(20) 



d9 



= 0, 



(I) 



= 0, 



(21) 



for antisymmetric modes. This is an eigenvalue problem 
and the spectrum tells us the stability of the solution. It 
is sufficient to examine the asymptotic behaviors of 6r 
for different values of A to extract sufficient information 



about the stability. Expanding Eq. ( 19 1 in small 9 gives 



d^Sr 
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IdSr 
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(22) 



which is independent of A and which possesses power-law 
solutions of the form Sr ^ 9^1"^, 9^1"^. These correspond 
to the symmetric and antisymmetric modes respectively. 

The asymptotic behavior in the opposite limit can 
be studied by making the transformations g(&) ~ 
6r{9)y/cot 9 and x = tan 6*, which results in 

d^g{x) sdg{x) 



dx^ 



+ p{x)^^+q{x)g{x) = 0, (23) 



dx 
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The asymptotic behaviors of these functions, as x 
+00, are 



(27) 
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and, 
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The asymptotic behavior of g{x) as x — > +00, for pos- 
itive values of A', can be computed by defining g(x) = 
exp{S{x)), where S{x) satisfies 



dx"^ 



dS 
dx 



dS 

+ p{x)—+q{x)=Q. (29) 



Using the eikonal approximation that S"{x) ^ (S"(x))^, 
which is valid for x +00, Eq. (29 1 can be solved 



asymptotically to give the two linearly independent solu- 
tions 
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where a series expansion in the form of. 
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(34) 



is performed to arrive at Eq. ( 33 1 



We see from the asymptotic formula, Eq. ( 33 1 , that 

1 



Sr2{x) = ^/xg2{x) ^1 + 



,5/2 



(35) 



asa:;-^ooor6'^7r/2. This means that Sr2{0) does 
not satisfy the boundary condition, Sr{9 = tt/2) = 0. 
The solution, 6ri{9), is the only solution that satisfies 
the boundary conditions, Eq. (20 1. 



To obtain the full eigenfunctions, we use the asymp- 
totic formula, Eq. (|32]) and ( 33 1 , as initial conditions 
and integrate numerically from a large value of a; = c 
(c = 10 in this case) back to a; = 0. The Gram-Schmidt 
orthonormalization procedure is employed to ensure the 
linear independence of the two eigenfunctions. The eigen- 
functions are normalized such that 



Sri{x)Srj{x)dx — Sij. 



(36) 
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FIG. 4: The eigenfunctions of Eq. ([T9| for A' = 0.1, 0.5, 1.0 
and 3.0. (a) The first eigenfunction satisfies the boundary 
conditions for symmetric modes, implying the instabihty of 
the dome solution, (b) The second eigenfunction does not 
satisfy the boundary condition at infinity. 



Fig. gshows Sri{e) and Sr2{0) for A' = 0.1, 0.5, 1.0 and 
3.0. From the graph, we confirm that ^1(6*) satisfies the 
boundary conditions, Eq.pO|), while ^2(6*) does not. 



Note that dri (6) satisfies only the boundary conditions 
for the symmetric modes, but not the anti-symmetric 
modes. We need a linear combination of Sri{6) and 
5r2{0) to form a solution that satisfies the latter. But 
since 6r2{0) does not satisfy the boundary condition at 
6 = 7r/2, such a linear combination would not satisfy it 
either. 



To conclude, there are always solutions to Eq. ( 19 1 sat- 
isfying the boundary conditions for the symmetric modes 
for every positive value of A, i.e., the domes are uncon- 
ditionally linearly unstable. This seems to be a contra- 
diction with the field observation of domes, which are 
presumably stable. We will postpone the discussion of 
this issue to the end of the next section, after we have 
discussed stalactite formation. 
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IV. STALACTITES 

In studying the formations of travertine domes near 
geothermal hot springs, it helps to study a similar geo- 
physical process, namely, the formation of stalactites, 
which are cylindrical structures formed by precipitation 
of calcium carbonate in limestone caves. Here, we will 
summarize the results Short et aZ. [51|7] obtained and ap- 
ply our formulation to study the stability of stalactites. 



A. Steady state 

As discussed earlier, the growth rate of stalactites is di- 
rectly proportional to the local fluid thickness, h. From 
the field observation, stalactite formation shares the fol- 
lowing features with dome formation: They both are cir- 
cularly symmetrical, formed under a shallow water lam- 
inar flow, and can be assumed to be locally flat. So, by 
using the analysis of dome formation, in particular, from 



Eq. ( 13 1, we have 



h = 



rsmi 



1/3 



(37) 



where f3 = 3i'Q/2ng is a constant. The dynamical equa- 
tion, Eq. ([2]), then becomes 



dn 



9^ 
902 



G 



r sin 9 



1/3- 



(38) 



where G depends on water chemistry and the input 
flux[6l |7j . Following the same strategy employed in the 
case of travertine domes, we obtain a uniformly translat- 
ing solution. 



r{9) 



ro 



sin 9 cos^ 9 ' 



(39) 



where the tip velocity vt comes in as an integration 
constant, and the scale rg = /3{G/vt)^. By defining 
p = r /tq, z = C/'f'a and using the trigonometric relation 
tan^ = —dz/dp, we obtain 



1 



= 0, 



(1 + z')^ P 
which is the result derived in Refs. OE]. 

B. Linear stability analysis 



(40) 



We study the stability of solution Eq. ( 39 1 by intro- 
ducing a perturbation: 



r{9) = f{9) + Sr{9)e 



At 



(41) 



where r is the unperturbed solution given by Eq. ( 39 ) 



and Sr is the perturbation. Substituting Eq. (41 1 into 



Eq. ( 38 I and expanding the resulting equation in 5r gives 
d^ T 



A' 



, dSr 
~d9 



cos ( 



1 



r sm f cos 



= 0, 



(42) 



where A' = ?>G^\/vf. We follow the same approach as in 
the case of the dome and study the asymptotic behaviors 
of the solutions of Eq. (42). For — > 0, we expand Eq. 



(42 1 in 6* and obtain 

, d6r 
~d9 ' 



A' 



1 



d^ 
dff^ 



96r = 0, 



(43) 



whose solution is given by r ~ where a = —1 — A. 
Because cr < for all A > 0, the solution diverges as 
9 Q. This shows that there are no eigenmodes for 
A > 0. As a result, we conclude that the steady-state 



solution Eq. ( 39 1 is linearly stable against the class of 



perturbations considered here. 

Let us also look at the asymptotics as a; 00 for 
completeness. Following the strategy employed in the 
study of dome stability, we make the transformation 
g{9) — tan0(5r(6') and x = tan^^. Eq. (42 1 then becomes 



where 



l9_ 
dx^ 



u(x) 



u(x) 



dx 



-8x 



v{x)g{x) = 0, 



A'(l + a;2)3/2 



and, 

v{x) 
As X ^ 00, 

u{x) ~ A'x^ -f 
and, 

v{x) 



A'(l + a;2)3/2 20x2 - 5 



(l + a;2)2 (i + a;2)5/2 



3A' 8 



3^ 

8x2 



O 



(44) 



(45) 



(46) 



(47) 



A'a 



2x 



20 



8x3 



O 



By following the same asymptotic analysis as we did in 
the last section, we get. 



gi{x) - exp 



-A'x3 3A'x 



and 



52(2;) 



1 

X 



10 

AV^ 



5A'x6 



O 



(49) 



(50) 



These can be used as the initial conditions to integrate 
numerically from a large value of x, giving the full eigen- 
functions. Again, the Gram-Schmidt orthonormalization 
procedure is employed. The two branches of solutions, 
'5?'i,2(6'), are plotted in Fig. [sl They do not satisfy the 
boundary conditions as they ooth diverge at = 0. So 
the stalactite solution is stable. 
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tan(e) 
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FIG. 5: The eigenfunctions of Eq. Q for A' = 0.1, 0.5, 
1.0 and 3.0. These solutions do not satisfy the boundary 
conditions, as they all diverges at 6 = 0. 



COMPARISON BETWEEN DOMES AND 
STALACTITES 



We have shown that there is a continuous spectrum 
of unstable modes for travertine domes, but stalactites, 
which are formed by an apparently similar process, are 
predicted to be linearly stable. We need to (a) explain 
why it is that domes can be observed in the field, and (b) 
interpret the source of the difference in stability between 
the two seemingly-related growth motifs. We initially 
found it surprising that there is a qualitative difference 
in stability, even though the dynamics of domes and sta- 
lactites seem to differ in only relatively minor ways: the 
growth of domes depends on the depth- averaged fluid ve- 
locity whereas the growth of stalactites depends on the 
fluid thickness. In both cases, the approximation of local 
flatness is used, so this is unlikely to be the source of the 
difference. 

Our interpretation is that the difference in stability 
arises from the direction of growth, and as a result, the 
manner in which surface tension effects correct the ze- 
roth order solutions we have discussed. The direction of 



growth is important, because it dictates the way in which 
shape perturbations propagate. For domes growing with 
sufficiently large Vt, shape perturbations are advected 
away from the vent down the body of the dome, in a man- 
ner reminiscent of the way in which shape perturbations 
are advected down the body of a growing dendrite |35j. 
These perturbations may also grow during this process, 
but the development of this instability is in practice reg- 
ularized by any non-zero surface tension, leading to con- 
tact line formation, film break-up and the formation of 
rivulets. This heuristic argument is supported by the 
shape of the linear stability eigenfunctions shown in Fig. 
|4] For stalactites, on the other hand, the fluid becomes 
increasingly thick as it flows down towards the tip, and 
perturbations only increase the growth velocity of the 
tip, rather than cause growing instabilities away from 
the tip. Thus, the only place where surface tension is 
significant is at the tip of the stalactite, where the sur- 
face tension holds a water droplet until the droplet be- 
comes too heavy and drops. This dynamics, we believe, 
mainly contributes to the precipitation rate at the tip, 
which affects only the growth rate of the whole stalac- 
tite. In other words, it only renormalizes the value of vt, 
which, in any case, is a fitting parameter. Surface tension 
is, therefore, not important in the dynamics of stalactite 
formation and it should not affect its stability. 

Returning now to the case of travertine domes, we con- 
clude that the unstable modes are small near the vent 
and grow in amplitude near the tail of the dome. This, 
however, is precisely the region where the film becomes 
thin and contact line formation can occur, leading to the 
fiuting pattern observed in the real systems. The pre- 
cipitation rate in this region is also lower, due to the 
depleted Ca^+ concentration, and this helps stabilizing 
the domes too. It is possible that the growth of the in- 
stabilities predicted here triggers the formation of con- 
tact lines and film break-up. Thus, we conclude that the 
dome is in some sense similar to the problem of dendritic 
growth, where a smooth tip is followed by a train of side- 
branches, widely interpreted to be due to a noise-induced 
instability [36, 37 . It is possible that the full inclusion of 
surface tension in the model would have as important 
a role in selection and stability as it does in dendritic 
growth [351 [35] ■ 



VI. DAMMING INSTABILITY 

Having studied the formation of domes and stalactites, 
we now try to understand some aspects of the large scale 
morphology of hot spring landscapes. We see in Fig. [TJa 
that the pattern formed is complicated, with ponds of 
similar shapes but different sizes. Empirical data shows 
that the distribution of pond sizes indeed follows a power 
law|20j. This scale-invariant pattern hints at an underly- 
ing scale-invariant precipitation dynamics, i.e., adynam- 
ics without a characteristic length scale. 

It is difficult to predict analytically the statistical prop- 
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erties of the landscape, such as the pond size distribution, 
due to the mathematical complexity of the equations in- 
volved. We can, nevertheless, study a simple case of pre- 
cipitation over a planar slope. By studying the linear sta- 
bility of this dynamics, we should be able to expose the 
essential physics of the formation of these scale-invariant 
patterns. The nonlinear regime of the modeling can be 
studied using the cellular model we introduced earlier. In 
this section, we consider a one-dimensional flow down an 
inclined plane, and evaluate the linear stability spectrum. 

The fluid flow in travertine systems is, unlike in the 
cases of dome and stalactite formations, generally turbu- 
lent. It is therefore necessary to use the formulation of 
Eq. ([7|-pT|). The turbulent drag leads to a steady flow 
regime, about which we linearize. Since the angle 9 is 
the same along a constant slope, it is more convenient to 
use the arc length, s, as the independent variable in the 
growth equation, so the dynamics of local curvature, k, 
is given by 



(51) 



where the subscript n denotes a derivative taken at a 
point moving along the outward normal of the curve. 
This, together with the Dressier equation, Eq. ([7|-pT|, 
gives the complete description of the system. 

We scale the independent variables to their natural 
units. 



9 = 



ln(l — (jnh), 



(58) 



where we deflned the Froude number, Fm = /gR. 
The uniform solution of this set of equations is given 

by 



Mo 



/ sin 61 

CfFra ' 



/l= 1, 



K= 0, 



(59) 



(60) 



(61) 



(62) 



where is the initial inclination of the slope. The lin- 
ear stability analysis is performed by adding harmonic 
perturbations to the solution. 



Wo 



Oo + ^uoe'^*+^*, 



s+\t 



(63) 



(64) 



(65) 



t' 



U 

—i 

R 



s 



c 

R' 



and dcflnc the following dimcnsionlcss variables. 



uo 



h' = 



k! = Rk, a — 



H 

R' 



(52) 



(53) 



where U, H and R arc the characteristic scales of the 
fluid velocity, fluid thickness and the landscape respec- 
tively, and a is the ratio between the H and R, which is 
small in the regime of shallow water flow. The governing 
equations then become (we drop all the primes on the 
variables for simplicity). 



dt 



Guo, 



(jF, 



dua dE _ ~CfF,nul 



dt ds h{l- ^) ' 



1 - (TAv/l a— + 7^=0, 

at OS 



with 



C + ahcos9 



Ph 



aFmul 



pg 2{l- anhy 



(54) 



(55) 



(56) 



(57) 



OS 



(66) 



and linearizing the resultant equations to the first order 
in the perturbations, resulting in three equations for Suq, 
6h and S9, 



ipXS9 = p^GSuq 



(A + ipuo)Sh + ipSuo 



69 = 0, 



(67) 



(68) 



aFmXSuo = S9{cos 9 + ipa sin 9 + p^a^u^Fm) 
^^ CfFmulaip 

+Sh{~ipa cos 9 + CjFmul) 
+Suo{-ip(TF„iUa ~ 2CfF„rUa) (69) 

A single dispersion relation can be obtained by com- 
bining all three equations and eliminating Suq, Sh and 
69. The result is a cubic equation in A, 

+ a2ip)X^ + ai{p)X + aoip) = 0, (70) 



where 



a2 (p) = 2iuop + 



2CfUo 



(71) 
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FIG. 6: The damming instability spectrum with parameters 
ieo,G,Fm,Cf,a) = (7r/6, 10"*, 10,0.1,0.01). (a)-(c) The real 
parts of the three branches of solutions. The first branch, 
Ai, is positive for all p, implying that the solution is uncon- 
ditionally linearly unstable, (d) The imaginary parts of the 
solutions. 



ai(p) 



' G sin ( 

^ m 

' iG COS ( 



+p 



CfulG 
2 



cos ( 

m. 



(tF„ 



(72) 



4 , _3 , auoGcos9 
ao(P) = P I -crwnG 



"■0 

-iGuQ 



sm 
-Guq cos 9 



2F„ 



0F^ + iGCful 



(73) 



VII. CONCLUSION 



By combining fluid dynamics and surface growth kine- 
matics, we formulated a mathematical framework to 
study geological pattern formation due to carbonate pre- 
cipitation and applied it to study the formation and sta- 
bility of a variety of motifs. The theory successfully pre- 
dicted the shape of observed spherically symmetric domes 
for angle 9 less than a critical angle 9c- By compar- 
ing with results from a cellular model, we showed that 
the departure of our theoretical prediction from observa- 
tion for 9 > 9c is due to the neglect of surface tension. 
We also showed that domes are linearly unstable towards 
axisymmetric perturbations, but the instability is mani- 
fested in the tail of the dome away from the vent. The 
instability is masked by the thinning of the fluid fllm 
and ultimately the formation of contact lines due to sur- 
face tension. This contrasted with the case of stalactites, 
whose growth forms are linearly stable to axisymmetric 
perturbations. The difference between the stabilities of 
the dome and stalactite solutions is attributed to the dif- 
ferent geometries and the different role surface tension 
plays in these two systems. 

This formulation cannot predict the complex landscape 
formed in the fully nonlinear regime, but a linear sta- 
bility analysis for a one-dimensional flow showed that 
the apparent scale-invariant landscape is consistent with 
our equations. In future work, we hope to examine the 
full two-dimensional instability problem, in order to in- 
vestigate the dynamics of pond formation, possibly as a 
transverse morphological instability, akin to meandering 
in step-flow processes on vicinal surfaces ^D]. 



For the parameter set {9o,G,Fm,Cf,a) — 
(7r/6, 10~^10,0.1,0.01), the three roots of the Eq. 
( |70| , Ai, are computed numerically and are plotted in 
Fig. |6] From the graph, we see the flrst branch of the 
solutions is always unstable, while the remaining two 
branches are always stable, implying that the solution is 
unconditionally linearly unstable. This is the damming 
instability. 

To conclude, we found that the trivial flow down a con- 
stant inclined plane is unstable towards all length scales, 
suggesting that when fully developed into the nonlinear 
regime, the landscape would have no selected length scale 
- a surmise in accord with field observations and our cell 
dynamical system simulations. 
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